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We develop a time-dependent theory to investigate electron dynamics and photoionization pro¬ 
cesses of diatomic molecules interacting with strong laser fields including electron-electron correlation 
effects. We combine the recently formulated time-dependent generalized-active-space configuration 
interaction theory [D. Hochstuhl and M. Bonitz, Phys. Rev. A 86, 053424 (2012); S. Bauch et al., 
Phys. Rev. A 90, 062508 (2014)] with a prolate spheroidal basis set including localized orbitals and 
continuum states to describe the bound electrons and the outgoing photoelectron. As an example, 
we study the strong-field ionization of the two-center four-electron lithium hydride molecule in differ¬ 
ent intensity regimes. By using single-cycle pulses, two orientations of the asymmetric heteronuclear 
molecule are investigated: Li—H, with the electrical field pointing from H to Li, and the opposite 
case of H—Li. The preferred orientation for ionization is determined and we find a transition from 
H—Li, for low intensity, to Li—H, for high intensity. The influence of electron correlations is studied 
at different levels of approximation, and we find a significant change in the preferred orientation. For 
certain intensity regimes, even an interchange of the preferred configuration is observed, relative to 
the uncorrelated simulations. Further insight is provided by detailed comparisons of photoelectron 
angular distributions with and without correlation effects taken into account. 

PACS numbers: 33.80.Eh,31.15.vn,31.15-p 


I. INTRODUCTION 

The rapid progress in experimentally observing and 
even controlling electron dynamics in atoms and 
molecules demands powerful theoretical approaches; see, 
e.g., EH for reviews on this subject. One of the most 
challenging and, therefore, interesting tasks is the accu¬ 
rate description and understanding of the ultrafast and 
complex behavior arising from the electron-electron inter¬ 
action. It can be expected that mean field, i.e., Hartree- 
Fock-type approaches are insufficient, and that electronic 
correlations are important. These are especially diffi¬ 
cult to treat in time-dependent theories with more than 
two active electrons, due to the complexity of the multi¬ 
electron wave function and even more so if the continuum 
is included for photoionization; see Q for an overview. 

Of particular interest in the context of strong field 
physics are molecular systems due to their much more 
complex dynamics and degrees of freedom owing to their 
geometrical structure. With the development of align¬ 
ment and even orientation techniques [5H8J measurements 
in the molecular-fixed frame of reference become accessi¬ 
ble, which allows for an investigation beyond orientation- 
averaged quantities. 

One question is the preferred direction of electron emis¬ 
sion with respect to the electrical field direction of a lin¬ 
early polarized laser and the influence of correlation ef¬ 
fects in strong-field excitation scenarios of heteronuclear 
molecules. In a first approximation, the tunnel ionization 
maps the highest-occupied molecular orbital (HOMO) to 
the continuum H, and by recollision strong-field ioniza¬ 
tion was even used to illustrate the HOMO experimen¬ 


tally 0. 

However, experimental evidence using CO 
molecules a E2| showed that this simplified one- 
electron picture needs to be adjusted, as effects such as 
inner-shell polarizations fldl] , Stark shifts and orbital 
distortions [141 ] have impact on the ionization dynamics. 
The question to what extent electronic correlations are 
important remains open. In Ref. fl5| this topic has 
been addressed within a one-dimensional model of the 
four-electron LiH molecule, and a shift of the preferred 
direction of emission is observed when electronic correla¬ 
tions are included. The immediate question of whether 
these effects are present in a full three-dimensional 
analysis shall be answered by this work and completed 
by angle-resolved investigations. 

All of these above-discussed issues call for a gen¬ 
eral, time-dependent theory including external (possi¬ 
bly strong) fields beyond a perturbative approach. The 
fundamental equation describing the physics of these 
quantum systems is the (non-relativistic) time-dependent 
Schrodinger equation (TDSE). However, its direct nu¬ 
merical solution, even by means of supercomputers, is 
limited to systems consisting of only one or two elec¬ 
trons, e.g. helium [TgI |20| | or molecular hydrogen [2~i1. It3|. 
Semi-analytical theories, such as the strong-field approx¬ 
imation and tunneling theories [ 23 M 28 I] provide physical 
insight but often draw on a simplified picture of the 
electron-electron interactions. 

In order to solve the time-dependent Schrodinger 
equation for more than two active electrons includ¬ 
ing the electrons’ interactions, approximate numer¬ 
ical techniques need to be employed. These in¬ 
clude the time-dependent configuration interaction sin- 
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gles (TD-CIS) method [29M3li| . multi-configuration time- 
dependent Hartree-Fock (MC-TDHF) [32M36| or its gen¬ 
eralizations time-dependent restricted or complete ac¬ 
tive space self-consistent-field (TD-RAS/CAS-SCF) (37i - 
|40| and the state-specific-expansion approach [41] (see 
also Ref. [ 3 ] for an overview). Further, time-dependent 
density-functional theory (TD-DFT) [42|, [43j and time- 
dependent close-coupling solutions of the TDSE by us¬ 
ing pseudo-potentials for the description of more than 
two electrons fill . |45| have been applied to photoioniza¬ 
tion of molecules. Especially the MCTDHF family suffers 
from complicated non-linear numerics, and its applicabil¬ 
ity to photoionization is not yet fully understood. TD- 
DFT and the pseudo-potential approaches, on the other 
hand, rely strongly on the chosen functionals or poten¬ 
tials with unknown accuracy and lack tunable parame¬ 
ters to achieve convergence to the fully correlated solu¬ 
tion. One of the most successful methods which bears 
some similarities to our present approach, is the time- 
dependent R-matrix method I f til I l8l |. 

The aims of the present work are (i) to provide a fully 
ab-initio time-dependent approach to electron dynam¬ 
ics in diatomic molecules exposed to strong laser fields 
including a systematic (i.e. controllable) approach to 
electron correlation without relying on pseudo-potentials 
and (ii) to demonstrate the method by shining light 
onto the question of whether electronic correlation de¬ 
cides from which end an electron leaves a heteronuclear 
molecule which is exposed to a strong electric single-cycle 
pulse. Our approach is based on the time-dependent 
generalized-active-space configuration interaction (TD- 
GAS-CI) formalism which we apply within a prolate 
spheroidal single-particle basis set in combination with 
the well-established partition-in-space concept to tackle 
the scattering part of the Hamiltonian. 

The paper is organized as follows. After a brief intro¬ 
duction into the theory of TD-GAS-CI, we give a detailed 
overview on the used basis set and details of our imple¬ 
mentation, in Section eh Technical aspects and the ex¬ 
plicit formulas and strategies of their efficient numerical 
handling are presented in the corresponding appendices. 
In sections Hill and ITVl we show illustrative numerical ex¬ 
amples and demonstrate the abilities of the present ap¬ 
proach. We focus on the LiH molecule in strong single¬ 
cycle infrared (IR) pulses and explore the influence of 
electronic correlations on the molecular photoelectron an¬ 
gular distributions (PADs) and the preferred direction of 
electron emission as a function of the geometrical set-up. 
The paper closes with conclusions and a discussion of 
future applications of the present theory. 


II. THEORY 

Let us consider N e \ electrons moving in the potential 
of two nuclei at positions Ra and Rb with charge num¬ 
bers Za and Zb- Throughout, we employ the Born- 
Oppenheimer approximation [49|, which decouples the 


nuclear and electronic degrees of freedom, and use atomic 
units ( m e = e = 47re 0 = 2|£ , R y( i| = a o = 1). The (elec¬ 
tronic) Hamiltonian is given by 


JVel JVel 

H(t) = YMt) ^ 


2=1 


i<j 


(1) 


with the one-electron contribution of the i-th electron 

hi(t) = ~V 2 i +V i + E(t)n, (2) 

consisting of the kinetic and potential energies with 


Vi = V{n) = 


Z A 

I ri - R a | 



(3) 


and the Coulombic electron-electron interaction. The 
time-dependent external laser field is denoted by E(t) 
and is included in dipole approximation using the length 
gauge via the position operator r t . 


A. TD-GAS-CI 

We solve the TDSE for the N e \ electrons, 

i^|VKt)) = ff(i)|*(i)>, (4) 

within the TD-GAS-CI framework 0, [l5|, 0,0| with the 
Hamiltonian 0. Thereby, we expand the many-particle 
wave function into a basis of time-independent Slater de¬ 
terminants 1$/), 

\m)= E ( 5 ) 

/£Vgas 

with time-dependent complex coefficients c/(t), which re¬ 
sults in the matrix representation of the TDSE, 

= E H u{t)cj{t) . (6) 

JGVqas 

The Slater determinants |$j) are constructed from single¬ 
particle spin orbitals x,(r, <T ) with the spatial coordinate 
r and the spin coordinate a and i — 1... 2 Nb, where Nj, 
is the dimension of the spatial basis set. Details of the 
orbitals are given in Sec. Ill Bl The matrix elements of 
the GAS Hamiltonian, Hu = can be evalu¬ 

ated either by directly using Slater-Condon rules (52j or 
by efficient techniques from (time-independent) quantum 
chemistry 0, 0. The most demanding task, besides 
the time propagation of Eq. m, remains the evaluation 
of the one- and two-electron integrals in Hij. Details 
of our strategy are given in Sec. Ill B 31 and App. EH The 
included determinants Vqas in sums (0 and 0 are cho¬ 
sen according to the GAS concept described in detail in 
Ref. [H}. Thereby, the method ranges from single-active 
electron (SAE) [sol , l54l . 0] to (exact) full CL 
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center A and center B and the electron, 
rA = \r- Ra/2\, 

r B = \r + R B /2\. (7) 

The vectors rq and Ri (i £ A, B) are given in Cartesian 
coordinates. The prolate spheroidal coordinates £ and 77 
are then defined as 


FIG. 1. (color online). Sketch of the coordinate systems for 
diatomic molecules [65] . The nuclei are labeled A and B, 
respectively. 


We propagate Eq. © using a short-iterative Arnoldi- 
Lanczos algorithm [561 . 157] which results in the repeated 
application of large-scale matrix-vector multiplications 
(up to 2 x 10 6 in one simulation) where the high degree 
of sparsity of the GAS Hamiltonian can be efficiently 
exploited. The algorithm is a ppl ied with an adaptive 
dimension of the Krylov space [56l . HIJ] for the propaga¬ 
tion of the wave packet with a time-dependent Hamilto¬ 
nian and with an adaptive time-step [59| for the propa¬ 
gation with a time-independent Hamiltonian after exci¬ 
tation with a laser pulse. 


e = !:d lr £ ’ (8) 

V=r - A 2^' ^[-1,1], (9) 

(j) = arctan ^ € [0, 27r]. (10) 

With this definition, center A is located at 2 < 0, and 
the binding potential for the electrons, Eq. ®, takes the 
form 

v = - * 2 A (Za + Z b )Z+(Z b -Za)v]- ( 11 ) 

a(r - V ) 

Explicit expressions needed for the implementation are 
comprised in the Appendix. 


B. Single-particle basis 


For the efficient solution of Eq. ©, a proper single¬ 
particle spin-orbital basis ( r , a ) with an associated spa¬ 
tial orbital basis <pi(r) is required to construct the de- 
terminantal basis |$j). For quantum-chemistry calcu¬ 
lations, typical basis sets are founded on expansions in 
localized functions, such as Gaussian- or Slater-type or¬ 
bitals (see, e.g., [53|)- These sets achieve a high precision 
for bound-state properties but lack an efficient descrip¬ 
tion of the scattering part of the Hamiltonian. For atomic 
targets, i.e., single-center potentials, typically mixed ba¬ 
sis sets with radial grids in combination with spherical 
harmonics are used 3113 . However, for multi-center 
geometries, the convergence exhibits unfavorable scaling 
properties for our case with the required expansion of 
the angular coordinates To overcome this prob¬ 

lem for diatomic molecules, we use two-center prolate 
spheroidal coordinates, in the following. Alternatively, 
also for larger molecules with more complicated geome¬ 
tries, a combined Gaussian and discrete variable repre¬ 
sentation (DVR) can be applied [62j-[64|]. 


1. Prolate-spheroidal coordinates 

The two-center problem of the diatomic system can 
be handled efficiently by prolate spheroidal (confocal el¬ 
liptic) coordinates. Therein, we define a = R/2 = 
\Ra — Rb\/% as half the distance between the two centers 
(nuclei), see Fig. [I] ta and r B are the distances between 


2. Spatial Basis 


For the set-up of the spatial basis, we follow closely 
Refs. [(31;, 661468] . We use a direct product basis where 
£ is represented by a Finite-Element DVR (FEDVR) ba¬ 
sis whose flexibility regarding the density of grid points 
avoids complicated coordinate scalings [69l . 17 O] . Coordi¬ 
nate 77 is handled by a usual Gauss-Legendre DVR [13, 
which is well suited for this problem because spheroidal 
wave functions are represented by Legendre polynomi¬ 
als. The spheroidal wave equation is similar to the one- 
particle Schrodinger equation (see, e.g., 0 ). 

Although L 2 (L is the electronic orbital angular mo¬ 
mentum) does not commute with the Hamiltonian, A 
(component of the electronic orbital angular momen¬ 
tum along the internuclear axis) does, and the associated 
quantum number m is a “good” one [72] , We, therefore, 
expand the (/^dependent part of T into the eigenvectors 
of A, which, in our case, outperformed a Fourier-Grid- 
Hamiltonian basis in the <fi coordinate [42j |: 




- 7 = (£, 77 ) exp (im^). (12) 

V Z7T 
v 


Note that in Refs. [61, [63 spherical harmonics 
Yi m [arccos(r 7 ), d>] for the 77 and (f> coordinates are used 
[73] . This basis shows slightly better convergence than 
a DVR in 77 . However, the resulting electron integrals 
are less sparse and the basis is non-orthogonal, which 
complicates the determinantal basis. 

To fulfill proper boundary conditions, we use Gauss- 
Radau quadrature for the first finite element, which en¬ 
sures that there is no grid point at the singularity £ = 1. 








4 


Gauss-Lobatto quadrature is used for the remaining ele¬ 
ments, as usual in FEDVR 0- To avoid singularity of 
the kinetic energy matrix, i.e., to render the matrix in¬ 
vertible for the calculation of the two-electron repulsion 
integrals (see App. [A}, the very last DVR point of the 
last element is not included in our grid. Thereby an in¬ 
finite potential barrier at the grid end is created, which 
forces the wave function to vanish asymptotically (Dirich- 
let boundary condition). If the grid is large enough (i.e., 
reflections are avoided), effects due to this procedure are 
negligible. 

The one-electron primitive functions used in this work 
are 


in Ref. [il|, the Hartree-Fock problem for the N e \ — 2 
electronic problem is solved with the exchange-potential 
included, which is appropriate for obtaining localized vir¬ 
tual orbitals. The outer part, £ > £ s , of the basis consists 
of non-rotated, “raw” functions which describe the wave 
packet in the continuum accurately. One block of the 
coefficient matrix C is hence diagonal, see Appendix lB 21 
An exploitation of the properties of this basis is inevitable 
for a fast and memory-friendly code [recall that the two- 
electron integrals scale, for an arbitrary basis, as 0 (N£) 
whereas for the DVR basis set they scale as 0{N%)\. Ap¬ 
pendix [B] gives details for the efficient transformations of 
the integrals. 


/»£(£,» lA) = 


a 3 (C 2 -^) 


«XT(^, (13) 


where a multi-index k has been defined for convenience. 
The form of the functions 9(x) and the matrix elements 
of the kinetic, potential and interaction energies are given 
in App. [A] along with details of their derivation. 


C. Observables 

In this part, we discuss the extraction of the relevant 
observables from the GAS wave function in the mixed- 
prolate basis set. More basis-independent details can also 
be found in Ref. fl5| . 


1. Angular distributions and photoionization yields 


3. Partially rotated basis 

In analogy to Refs. [l5|, [5(1, we use a partition-in- 
space concept to allow for an efficient description of the 
photoionization process. Similar strategies are also ap¬ 
plied in time-dependent 17-matrix theory, e.g. [47], and 
in Ref. (74j . Here the basis set is split at £ = £ s into 
two parts: an inner region, £ < and an outer part, 
£ > £ s . The splitting point £ s is chosen such that it 
coincides with an element boundary of the FEDVR ex¬ 
pansion. This assures the continuity of the wave function 
across the grid and avoids the evaluation of connection 
conditions, see Ref. 0 for a detailed investigation in one 
spatial dimension. 

The basis in the inner region, £ < £ s , is constructed 
from Hartree-Fock-like rotated orbitals, 


Angular distributions of photoelectrons contain a 
wealth of information (see, e.g., 0 ), and, especially 
in strong fields, dynamical properties of the rescatter¬ 
ing process lead to rich structures (7(|. The molecu¬ 
lar angular-resolved photoionization yield (photoelectron 
angular distribution, PAD) is defined as 

/>27T /»oo 

y(8,t) = / / dcj)dr r 2 p(r, 9, <(>; t). (15) 

Jo Jr c 

p is the charge density in spherical coordinates (9 is the 
azimuthal angle in the 2 -a:-plane), 

N b 

P(r;t) = J2 D ki( t )fk(r)*fi{r), (16) 

kl 


= ^ C ijfA r ), ( 14 ) 

3 

where C is the orbital coefficient matrix. This rota¬ 
tion is needed because the energy of a truncated Cl 
wave function changes under a unitary transformation 
of the underlying single-particle basis; a good single¬ 
particle basis drastically enhances convergence with re¬ 
spect to the size of the truncated Cl space, but comes 
at the cost of expensive integral transformations which 
destroy the desired (partial) diagonality of the integral 
matrices [TB!, [53]. Since we are interested in one-electron 
photoionization with the simultaneous excitation of the 
ion, we follow the detailed investigations in [Ie| and use 
pseudo-orbitals based on the A^ e i — 2 electron Hartree- 
Fock problem for the virtual orbitals in the rotated part 
of the basis. Here, in contrast to the procedure shown 


with the spin-summed single-particle density matrix 

= (17) 

a 

a' ka . and dka are the creation and annihilation operators, 
respectively, of a spin orbital with spatial index k and 
spin a. The critical radius r c is chosen to be sufficiently 
large such that only the “ionized” part of the charge den¬ 
sity is used for integration. This neglects the long-range 
character of the Coulomb potential and is strictly valid 
only for r c —> oo. Therefore, several r c are used to check 
convergence. r c lies usually in the outer region of the 
partially rotated basis. 

The orbitals sampled at grid points in spherical coor¬ 
dinates, fk{r,9,(f >), cf. Eq. (1161) . can be efficiently stored 
in a sparse vector format by exploiting the locality of the 
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FEDVR functions. This decreases the memory require¬ 
ments and the computation of the charge density (1161) in 
typical computations by more than three orders of magni¬ 
tude. The photoionization yield at time t can be retrieved 
from the integrated PAD: 

V{t)= f d0y(9, t) sin 9. (18) 

Jo 


2. Photoelectron energy distributions 

The momentum distribution of the photoelectron is ob¬ 
tained by using the Fourier-transformed basis functions 
[H]. Only basis functions outside r c are used, ignoring 
the centrahregion. This approach is exact for sufficiently 
large r c |77fl. The Fourier transform of function ft. is 
defined asll| 

fkip ) = J dr exp (—ip • r)f k (r), (19) 

where r and p are vectors in prolate spheroidal coordi¬ 
nates and p ■ r is the inner product in these coordinates. 
If the ^-component of p is either 0 or 7r, analytical expres¬ 
sions of Eq. (fTTJl) can be obtained. However, the integral 
kernel is nonanalytic which prohibits the usage of the 
DVR properties (Gauss quadrature) of the basis func¬ 
tions for the integration. Therefore, we employ the Fast 
Fourier Transform of the basis functions in Cartesian co¬ 
ordinates. An application of Eq. m with the Fourier- 
transformed basis functions gives then the momentum 
distribution. 


III. APPLICATION TO LITHIUM HYDRIDE 



FIG. 2. (color online). Isosurfaces of the restricted Hartree- 
Fock HOMO (orbital to the right, green) and the core orbital 
(blue, to the left) of the LiH molecule. The Li nucleus is 
located to the left (this corresponds to configuration Li—H in 
Sec. [JV]). The iboview program was used for the generation 
of the orbitals using an isosurface with a density-threshold of 
90.99% M- 

and of sin 2 shape, 

^W={ 0 E ° Sil,(f ) 2cOSM) ’ (2D 

with the photon energy uj , the amplitude Eo, and the 
carrier-envelope phase (CEP), <^cep- 

The simulations are carried out within the fixed-nuclei 
approximation, which is well-justified since the dynamics 
of the nuclei are on a much longer timescale than the 
considered pulse durations. All data are retrieved by 
using the length gauge, cf. Eq. ©, which is preferable 
over the velocity gauge in the case of tunnel-ionization 
dynamics with few-cycle pulses (83j . 


Let us consider the diatomic molecule LiH, i.e., Za = 3 
and Zb = 1- It is the smallest ( N e i = 4) possible hetero- 
nuclear molecule and exhibits a spatial asymmetry with 
respect to its geometrical center. It is a frequently chosen 
theoretical model to test correlation methods (e.g., in one 
spatial dimension in Refs. mmm zsUsS). 

We use an internuclear distance of R = 3.015 (1.60 A), 
which is the equilibrium geometry at CCSD(T)/cc- 
pCV5Z level {§2] (without the frozen core approximation) 
and the experimental value [72|- The HF electronic struc¬ 
ture consists of a valence (HOMO) and a core orbital (Is 
of lithium), which are given in Fig. [21 

The electric field is linearly polarized parallel to the 
internuclear axis. We consider envelopes of Gaussian 
shape, 


Eg auss(^) — Eq exp 


r ( t-t 0 ) 2 1 

2a 2 


cos [u(t - to) + v?cep] , 


( 20 ) 


A. GAS partitions 

The TD-GAS-CI method is well suited for photoioniza¬ 
tion problems (l5l[50j . as it can be tailored to the problem 
at hand. For constructing the GAS, we assume that mul¬ 
tiple ionization is negligible due to the much larger ioniza¬ 
tion potential of LiH + . We consider three types of GAS 
partitions which are sketched in__Fig. [3J the single-active 
electron (SAE) approximation [sol |54] . [ssl] , the Cl singles 
(CIS) approximation and a complete-active-space (CAS) 
with single excitations from this subspace to the remain¬ 
ing orbitals in the outer region. Following Ref. [45j], we 
denote this type of GAS with C AS* with 

describing the number of electrons and v the number of 
spatial orbitals in the CAS. The star indicates the sin¬ 
gle excitations out of the CAS. For the LiH molecule, 
CAS*(2,p) describes a CAS with frozen core, and for 
CAS* (4,(2), all electrons are active. As a sidemark, we 
mention that this type of GAS is equivalent to a multi¬ 
reference CIS description. The size of the active space, 





6 


E 




im 


1 <M 

4^ 

I0O> 



GAS3 


} 


GAS2 

(SAE) 

GAS1 

(frozen) 


-f+ 

-f-+ 


GAS2 


GAS1 

(CIS) 


SAE 


CIS 


im 


{i^>} 
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FIG. 3. (color online). GAS divisions used in this work. Orbitals (labeled by rf>i) of different spin are assumed to have the same 
energy. The (red) arrows with dots show the allowed excitations. Striked-out arrows mean that no excitations are allowed. The 
nomenclature is SAE: single-active electron, CIS: configuration interaction singles, CAS: complete-active space with v spatial 
orbitals, the star indicates single excitations out of the CAS. 


i.e., v in CAS*(*,i/), remains to be chosen adequately 
and to be checked carefully for convergence. 

For LiH, the accurate description of correlation ef¬ 
fects requires orbitals with higher quantum numbers A 
in the <f> coordinate, i.e., m max > 0. In contrast to time- 
independent CAS-SCF simulations, where typically a few 
additional orbitals, such that open shells are filled in 
the CAS, give satisfactory results for the Stark-shifted 
ground-state energies, time-dependent calculations are 
more involved. Here, in order to describe intermediate 
states accurately, a larger size of the CAS is crucial. Rea¬ 
sonable CAS configurations contain closed subshells in 
the number of active orbitals, v, i.e., all orbitals with 
a certain symmetry. This leads, for example, to a rea¬ 
sonable space of CAS* (2,5) for orbitals with A = E and 
n-orbitals for two active electrons [CAS* (4,6) for four 
active electrons]. Typically, we use CAS spaces with or¬ 
bitals of up ton (m max = 1), CAS*(2,8) and CAS*(2,12), 
and A (m max = 2) symmetry, CAS* (2,10). 

The convergence of the method is illustrated with 
an example of the photoelectron angular distributions 
(PADs) for the strong-field ionization of LiH using single¬ 
cycle pulses in Fig. [I] see Sec. IIVI for the field param¬ 
eters. All GAS approximations predict the dominant 
electron emission in direction of the field polarization, 
0 = 0°, 180°. However, the SAE and CIS approximations 
drastically underestimate the total yield and fail to pre¬ 
dict the correct positions of the side maxima (see inset of 
Fig. lU. The strongest difference is observed in direction 
of 6 = 0°. By increasing the active space, a successive 
convergence is achieved and for a CAS* (2,8) only small 
differences appear in comparison to larger spaces (blue 
dashed vs. solid blue line). Therefore, we typically use 
the CAS* (2,8) model in the following. The convergence 
was checked for different field parameters additionally. 



FIG. 4. (color online). Convergence of the photoelectron 
angular distribution for different GAS at a field intensity of 
2.2 x 10 13 W cm -2 with m ma x = 1. Data is shown for single¬ 
cycle excitation of configuration H—Li, see Sec. lIVI for details. 
The PAD for CAS*(2,10) with m max = 2 (CAS*(4,4)) coin¬ 
cides with that of CAS* (2,8) (CAS* (2,3)) and is therefore not 
shown. 


B. One-photon ionization 

To demonstrate the method, we first consider the case 
of one-photon absorption in LiH. The parameters for the 
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mean photoelectron energy (a.u.) 

FIG. 5. (color online). Photoelectron energy spectrum for 
the ionization of LiH with a short (2.01 fs) pulse with 40.8 eV 
and an intensity of 0.088 x 10 13 W cm -2 within SAE and 
CAS* (2,8) (converged) approximation. 


electric field, Eq. ([20)1 . are w = 1.5 (40.8eV or 30 nm), 
E 0 = 0.005 (0.088 x 10 13 Wcm- 2 ), cr = 70 [2.82 fs full 
width at half maximum (FWHM) of intensity], to = 350 
(8.47fs), and y>cEP = 0. The TDSE in GAS approx¬ 
imation is propagated until f = 900 (22fs). We use 
up to three functions in (f> (m max = 1), ten in rj , and 
1386 in £. The inner region in £ consists of two elements 
with 10 and 18 basis functions each and ranges [1, 2] and 
[2,15), respectively. The non-rotated basis is formed by 
80 equidistantly distributed elements with 18 basis func¬ 
tions in each element within £ 0 uter £ [15, 800). 

The kinetic energy spectrum of the photoelectron is 
shown in Fig. [5] for the SAE approximation and the con¬ 
verged CAS*(2,8). The spectrum shows a strong peak 
at the expected position of E^ j n = w — Ip ^, where 
Ip ^ = |E’„| = 0.295 is the ionization potential accord¬ 
ing to Koopman’s theorem, i.e., the negative HF en¬ 
ergy of the valence orbital, E v . This peak is rarely 
shifted by correlations (solid line) but a series of addi¬ 
tional peaks appears at lower kinetic energies (see in¬ 
set of Fig. [5] E < 0.8). These can be attributed to a 
correlation-induced sharing of the photon’s energy be¬ 
tween the photoelectron and a second electron still bound 
in the ion (“shake-up” process). Thereby, the photoelec¬ 
tron energy is reduced and the ion remains in an excited 
state. The origin of this process is purely correlation- 
induced and can neither be described within neither the 
SAE approach, the CIS approach ca , nor TD-HF simu¬ 
lations [84J. The population dynamics of these states can 
be measured, e.g., by strong-field tunneling [85| . 
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FIG. 6. (color online). Photoelectron momentum distribution 
of LiH after one-photon excitation for a CAS* (2,8). An in¬ 
tegration cutoff of r c = 150 ao (radial coordinates) was used. 
See Fig. [5] and text for parameters. 


The corresponding angle-resolved momentum spec¬ 
trum depicted in Fig. (G) contains additional information. 
The photoelectron shows characteristic angular distribu¬ 
tions for the different peaks with distinct locations of the 
maxima. The outer circle with a radius of about 1.5 cor¬ 
responds to the main photoelectron peak at an energy 
around 1.2 in Fig. [5] The inner (fainter) circles stem 
from the shake-up state population and exhibit a signifi¬ 
cantly different angular dependence than the main peak. 
This is caused by the different selection rules for the si¬ 
multaneous excitation of two electrons and, therefore, 
the changed angular momentum of the escaping electron 
in comparison to the dominant ionization channel with 
LiH + in its ground state. 


IV. STRONG-FIELD IONIZATION OF 
LITHIUM HYDRIDE 


We now turn our attention to the case of strong and 
short pulses, for which a time-dependent theory is indis¬ 
pensable. Let us consider single-cycle pulses of the form 
of Eq. (1211) with uj = 0.057 (800 nm) which corresponds 
to a duration of 5.3 fs. 

The electrical field exhibits a strong CEP dependence 
with a predominating orientation at the maximum field 
strength, see Fig. [7] (b). This dependence corresponds 
to two different ionization scenarios, i.e., orientations of 
the molecule with respect to the field at the maximum 
intensity of the linearly polarized pulse. We will refer to 
these situations as Li—H, if the field points from the H 
to the Li end [panel (a) in Fig. [7] and H—Li, if the field 























FIG. 7. (color online). Sketch of the two ionization scenarios 
of LiH: Figure (a) and (c) show the potential at the peak of 
the field at to for the configuration Li—H (H—Li). The black 
horizontal line denotes the energy of the active orbital. Panel 
(b) shows the time dependence of the single-cycle pulse. 


points from the Li to the H end [panel (c)[. 

The single-particle basis is similar to Sec. IIII Bl but 
with 16 functions in ry and 920 functions in £: 8 and 
14 functions in the inner region and 100 elements with 
10 functions each in the outer region within £ 0 uter G 
[15,1000). The total number of basis functions in our 
simulation is 14720 (m max = 0), 44160 (TO max = 1) 
and 73600 (m max = 2). We further note by compar¬ 
ing the ionization potentials of LiH (Ip 1 ' 1 = 0.295) and 
LiH + (Ip 2 ' 1 = 0.825) that the Keldysh parameter [23| 
7 = \JU V I2I V is 1.7 times larger for the ion. Therefore, 
we conclude that single excitations into the outer region 
is a valid GAS approximation and double ionization is 
negligible. Convergence aspects of the size of the CAS 
were discussed in Sec. eeh 

The charge densities at different times during the 
simulation for a field strength of Ei = 0.025 (2.20 x 
10 13 W cm -2 ) and the two scenarios Li—H and H—Li are 
given in Fig.[5]for the converged CAS*(2,8) (rn max = 1). 
The main dynamics happens after about 3 fs when the 
maximum amplitude of the pulse is reached and the tun¬ 
neling ionization sets in. The ejected part of the wave 
packet exhibits a characteristic angular distribution, vis¬ 
ible in the logarithmic density plot. After the ionization, 
the molecular ion remains in Rydberg states and per¬ 
forms coherent oscillations between the electronic states. 
A closer inspection of the released wave packet after the 
pulse is over (t = 13.5 fs) reveals a significant difference 
for scenarios Li—H and H—Li in both, the angular distri¬ 
bution and the absolute yield. The results presented in 
Fig. [5] confirm also the observation in Ref. [15| utilizing 
a one-dimensional model for LiH that ionization for field 
strength Ei is preferred from the Li-end, i.e., using the 
configuration Li—H. This finding will be quantified and 
discussed in detail in the remaining part of the paper. 


A. Orientation dependence of electron emission 

To address the question of whether the ionization yield 
is larger if the linearly polarized light is pointing from 
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FIG. 8. (color online). Snapshots of the charge density (log¬ 
arithmic plot) for the two orientations of the LiH molecule 
at different times during the excitation with the single-cycle 
pulse using an intensity of 2.20 x 10 13 W cm' 2 for both sce¬ 
narios, Li—H (left) and H—Li (right), cf. Fig. [3 The classical 
force F = — W points downwards. 


the Li- to the H-end (configuration H—Li) or vice versa 
(Li—H), the ratio 


C = lim 


V 


H-Li 


(*) 


t —^oo pLi ^(t) 


( 22 ) 


is defined [Ij|. ? is smaller (larger) than one if the elec¬ 
tron is ejected mostly in the direction from the H- to 
the Li-end (Li- to H-end); see also Fig. [7] In Ref. [l3| a 
similar parameter was defined for the CO molecule and 
a change of the direction for single-active-orbital calcu¬ 
lations in comparison to (uncorrelated) TD-HF calcula¬ 
tions including inner orbitals was found. 

Let us first consider a field strength of E\ = 0.025 
(2.20 x 10 13 Wcm -2 ). From the previous discussion, we 
expect c < 1, since ionization is preferred for configura¬ 
tion Li—H, cf. Fig. [5] This can be understood due to 
the increased total density at the Li nucleus and the di¬ 
rection of the classical force F = —X/V acting on the 
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Method 

^max 

c 

SAE 

0 

0.30 

CIS all electrons active 

0 

0.32 

CIS frozen core electrons 

0 

0.31 

CAS*(2,3) 

1 

0.31 

CAS*(2,5) 

1 

0.35 

CAS*(2,8) 

1 

0.43 

CAS*(2,12) 

1 

0.44 

CAS*(2,10) 

2 

0.44 

CAS* (4,4) 

1 

0.31 


TABLE I. Ionization asymmetry parameter Eq. (1221) . 
for the single cycle pulse with field strength E\ (2.20 x 
10 13 Wcm -2 ) calculated in different GAS approximations 
with one [SAE/CIS], two [CAS*(2,») ] and all four [CAS*(4,») 
] electrons active, c < 1 corresponds to a preferred configura¬ 
tion Li—H for ionization. 


electronic density. For the opposite configuration, H—Li, 
large portions of the electronic density have to “pass” an 
additional potential well, cf. Fig. [3 This simple picture 
is confirmed by SAE and CIS calculations in Tab. Q] with 
values of c < 0.3. The total ionization yield is, therefore, 
about a factor of three larger for Li—H than for H—Li at 
field strength E\. 

Additionally, Tab. [Q demonstrates the behavior of <; 
with respect to electron correlations by using different 
sizes of the CAS. Similar convergence behavior is also 
found for other field strengths. By successively increas¬ 
ing the CAS, first by including only two active elec¬ 
trons, CAS*(2,*), c increases, as it was observed in one¬ 
dimensional LiH [lH]. Most of the correlation contri¬ 
butions are captured by a CAS*(2,8) with a value of 
C = 0.43 and only less than 4% change is found by further 
increasing the active space. SAE and the CIS approxi¬ 
mations, however, predict values between 0.30 and 0.32. 
Therefore, correlations shift the preferred end of ioniza¬ 
tion from the Li to the H end significantly. 

We note that too small CAS, e.g. CAS*(2,3), give in¬ 
accurate results because of a bias due to an improper 
selection of additional important configurations, which is 
a general pitfall of multi-reference methods |8g|. How¬ 
ever, to test the frozen-core approximation (correla¬ 
tions arising from the two core electrons are not taken 
into account), we also performed calculations with all 
four electrons active for a small active space [CIS and 
CAS*(4,4) in Tab. [I]. By comparing to CIS (frozen core) 
or CAS*(2,3), respectively, we find that the observable 
does not change substantially. For larger CAS, a simi¬ 
lar behavior is expected, as was demonstrated for one¬ 
dimensional systems in Ref. [l5|. Also an increase of 
m max from one to two does not change the result. For 
both, the SAE approximation and CIS, even ro max = 0 
is sufficient, because the used orbitals in <f> are the corre¬ 
sponding eigenvectors of the one-electron problem. 



FIG. 9. (color online). Ratio c, cf. Eq. (1221) . of the pho¬ 
toionization yields for the two configurations [see Fig. [7] as a 
function of the field intensity for the relevant GAS approxi¬ 
mations. ? < 1 (? > 1) corresponds to ionization from the Li- 
(H)-end. See text for all other parameters. 


B. Intensity dependence 

In the previous paragraph, we found a dominating re¬ 
lease for configuration Li—H at a rather high intensity 
of 2.20 x 10 13 Wcm~ 2 , which is well above the barrier. 
This is understandable by the higher electron density at 
the Li nucleus. The HOMO, however, is more localized 
at the H end (green orbital to the right in Fig. [3, in 
contrast to the full single-particle density. Thus, enter¬ 
ing the tunneling regime by reducing the intensity, the 
mapping mechanism of the HOMO to the continuum be¬ 
comes important. Therefore, we expect larger values of c 
for low intensities, and thus by increasing the intensity, a 
decrease of <;; in other words, we expect a shift from the 
preferred H—Li to the Li—H configuration. 

This expectation is readily verified by the simple SAE 
approximation in Fig.[9](red line with squares). For small 
intensities, ? > 1, exhibiting values of about 1.06 indi¬ 
cating a slightly more favorable ionization of the H—Li 
configuration. For this situation, the CIS approximation 
(bright yellow line with triangles) fails to describe the 
change of preferred configuration. A c of maximal 0.6 is 
predicted by CIS. The better qualitative description of 
the physics by the SAE approximation compared to CIS 
is probably due to error cancellation. 

By increasing the intensity, a monotonic transition 
to the previously discussed strong-field over-the barrier 
regime with ? < 1 is observed. The change of the pre¬ 
ferred configuration for ionization from H—Li to Li—H 
(c « 1) occurs for the SAE approximation at around 
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1 x 10 12 Wcm -2 . For larger intensities, SAE and CIS 
calculations approximately coincide and the curves ap¬ 
proach a ratio of about 0.3, i.e., a factor of three higher 
yield for Li—H. For very high field strengths, the value 
increases again, which results from the fact that for ex¬ 
tremely intense pulses, the asymmetry of the binding po¬ 
tential is negligible in comparison to the excitation poten¬ 
tial and, therefore, a ratio of ? = 1 in the limit Eq —> oo 
is expectable. 

Electronic correlations [CAS*(2,8), blue line with cir¬ 
cles] change this picture qualitatively, in particular for 
smaller intensities: At a specific intensity range from 
around 3 x 10 12 Wcm -2 to 1 x 10 12 Wcm -2 , the corre¬ 
lation contributions interchange the dominant direction 
of emission, <; > 1 for CAS* (2,8) and <; < 1 for SAE 
and CIS. At 1 x 10 12 Wcm -2 , <, is much larger (1.64) 
for CAS*(2,8) compared to SAE (1.06), indicating that, 
during the slow tunneling process of the electron from 
the H-end to the continuum, much electron correlation 
can be built up. For large intensities, a similar trend as 
for the uncorrelated calculation is found, but with larger 
absolute values of ?, similar to the case of field strength 
Ei , which is discussed in detail above; see Tab. [TJ 


C. Photoelectron angular distributions 

We now turn our attention to the fully resolved PADs 
which contain more information than the integral quan¬ 
tity <r. The correlated CAS* (2,8) PADs for the two molec¬ 
ular orientations with respect to the electric field at dif¬ 
ferent field strengths are depicted in Fig.flOlin polar plots. 
As expected, the dominant ionization occurs along the 
field polarization axis for all considered intensity regimes. 
However, for all field intensities shown, the most favor¬ 
able ionization direction is the opposite direction of the 
field, see sketch at lower right of Fig. [TIB Further, the 
ejection direction of electrons differs from the preferred 
configuration for ionization measured by the ratio of the 
integral quantity ?. For example, at the highest field 
strength (top left in Fig. [TUI) the PAD shows a prefer¬ 
ence of H—Li over Li—H, whereas the value of ? < 1 
prefers Li—H [see Eq. (EU1) and Fig. [9j. This can be at¬ 
tributed to the propagation of the released electrons in 
the field: whereas ? depends mainly on the orientation 
of the molecule with respect to the main peak of the 
field, the PADs are modified by all cycles of the pulse 
and the cycles following the field maximum at which 
the tunneling-release of the electrons predominantly oc¬ 
cur can change the direction of electron emission signifi¬ 
cantly. This is similar to the rescattering mechanism for 
higher-harmonics generation and above-threshold ioniza¬ 
tion. This picture is verified by the time-dependent elec¬ 
tron densities in Fig. [5] where after the main peak of the 
pulse at t = 3.4 fs, the electron density gets accelerated 
by the smaller side-extremum in opposite direction (com¬ 
pare with t = 4.4fs in Fig. [Hj). A similar effect is, e.g., 
observed in strong-field ionization of atoms, where rescat¬ 


tering effects can drastically modify the angular distribu¬ 
tions of the photoelectrons (76|. 

With decreasing field strength (top left to bottom 
right), the shape of the curves along the main maxi¬ 
mum become more oblate and the smaller maxima first 
increase and then decrease again. For the highest field 
strengths, the PADs for the configuration H—Li (bright 
yellow curves) show also ionization contributions to the 
opposite direction pointing to the H end, which is not the 
case for Li—H (black curves). However, with decreasing 
field strength, the maximum pointing to the other di¬ 
rection is growing for configuration Li—H, and gets even 
larger than that of H—Li. Remarkably, the positions of 
the secondary maxima at this intensity are the same re¬ 
gardless of the position of the nuclei, but at intermedi¬ 
ate field strengths, an additional maximum for the H—Li 
configuration at the site of the H nucleus is visible. This 
indicates a higher angular momentum for the ejected elec¬ 
tron. 

To single out the influence of electron-electron cor¬ 
relations, the PADs from correlated CAS*(2,8) calcu¬ 
lation and those from a SAE calculation are shown in 
Fig. [TTJ For the Li— H configuration (left two columns 
in Fig. [TT1) . SAE (bright curve) underestimates the size 
of the side maxima, especially at intermediate field 
strengths (central panels). At smaller field strengths 
(0.351 x 10 13 Wcm -2 and 0.197x 10 13 W cm -2 ), the max¬ 
imum pointing away from the Li nucleus is either over- 
or underestimated, showing that no general pattern can 
be reasoned from correlation effects in PADs at different 
intensities. 

For the H—Li configuration and for all field intensities 
(right two columns in Fig. HD, SAE drastically under¬ 
estimates the side maxima at the site of the H nucleus. 
On the other hand, at intermediate intensities, the other 
side maxima are overestimated by the SAE approxima¬ 
tion. Thus, the favored direction of emission of electrons 
is decided by the electron-electron correlation for inter¬ 
mediate field strengths (see Fig. 0). 


V. CONCLUSIONS 

In this paper we presented a time-dependent approach 
to correlated electron dynamics following the excita¬ 
tion of diatomic molecules with strong electromagnetic 
fields. The method is based on the TD-GAS-CI approach 
using a prolate-spheroidal representation of the single¬ 
particle orbitals within a partition-in-space concept to 
allow for good convergence of the truncated Cl expan¬ 
sion. Thereby, parts of the multi-particle wave function 
close-by the nuclei are represented within a Hartree-Fock- 
like orbital basis and the ejected part is represented in a 
grid-like FE-DVR basis set. 

We illustrated the method by its application to the 
calculation of angle-resolved photoelectron spectra of 
the four-electron heteronuclear LiH molecule with and 
without taking electron-electron correlation contribu- 
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FIG. 10. (color online). Comparison of the PADs for the two configurations Li—H and H—Li at different field intensities for 
a correlated simulation with CAS* (2,8) using single-cycle pulses. All distributions have been scaled to fit in the range [0,1]. 
Note that the schematic (bottom right) gives the potential at the maximum of the electrical field. The dominant direction of 
the PAD includes also the propagation of the electrons in the field after ionization; see the text for a discussion. 
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FIG. 11. (color online). Molecular PADs at different field strengths showing the converged correlated simulations with a 
CAS*(2,8) and with the SAE approximation for the two orientations of the molecule (left two columns Li—H, right columns 
H—Li). The arrows indicate the predominating direction of electron emission. 


tions into account. To demonstrate the capabilities of the 
present approach, we then concentrated on the strong- 
field ionization of LiH using single-cycle pulses. The 
ionization yield for the two opposite orientations of the 
molecule along the linearly polarized electric field was 
calculated and an intensity-dependent shift of the pre¬ 
ferred configuration was observed: While for low intensi¬ 
ties in the tunneling regime, ionization for H—Li is larger, 


for high intensities well above the barrier, Li—H shows 
higher yields. In between both regimes, a smooth tran¬ 
sition is found. By turning on electronic correlations 
in the simulation, we find that especially yields in the 
tunneling regime are affected whereas the high-intensity 
regime is well described using the SAE or CIS approx¬ 
imations. Correlations shift the preferred configuration 
from Li—H to H—Li for low intensities and vice versa 
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for high intensities. In a certain intermediate intensity 
regime even an interchange of the preferred configura¬ 
tion in comparison to uncorrelated calculations is ob¬ 
served. Additionally, angle-resolved photoionization dis¬ 
tributions were presented and discussed, and the correla¬ 
tion effects were singled out by comparison to SAE cal¬ 
culations. 

Our results demonstrate the importance of electron- 
electron correlations in strong-field excitation scenarios 
of diatomic molecules. We expect the TD-GAS-CI ap¬ 
proach in combination with the prolate spheroidal ba¬ 
sis set to be applicable to larger systems such as the 
CO molecule and to arbitrary polarization of the excit¬ 
ing pulse in the near future, where experimental data is 
available 0,111. Further, the application to two-color 
excitation scenarios, such as streaking and XUV-XUV 
pump-probe, and the exploration of correlation effects in 
molecular systems on ultrashort time scales, e.g., post¬ 
collision interaction effects 0 . Isslj or the time-delay in 
photoemission is within reach. 
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the singularities at £ = 1, r) = ±1. Still, for odd to values, 
the exact eigenfunctions show non-polynomial behavior 
at these points and contain factors of (£ 2 — 1) 2(1 — 
r, 2 )^ 0,[H. Non-polynomial functions are poorly rep¬ 
resented by DVR, in which Gauss quadrature is used. 
Ther efore, for odd to, w e multiply the basis functions 9 n 

by V(£ 2 - !)/(£« - !) or \/(l - 77 2 )/(l ~ vl) t0 avoid a 

non-polynomial integrand [66] : 


N 


v i^n 

- Xi ’ 


(i, 


m even 

- y n (0 x | y 

€ 2 -i 

IJ-T’ 

m odd, 

(i, 


m even 

= Vn{v) X < [ 

T-r f 

m odd. 


(A4) 

(AS) 

(A6) 


{Vn) is the DVR grid point of the corresponding ba¬ 
sis function. These definitions differ in case of a bridge 
function in the FEDVR basis IIIlZQj]. 


1. Potential energy and dipole operator 


The potential, Eq. m, is independent of (j>. Hence, 
the potential matrix elements are evaluated using the 
DVR properties: 


Appendix A: Matrix elements 


Because, to our knowledge, the formulas of all needed 
matrix elements have not been stated in one single pub¬ 
lication or exhibit some misprints, we briefly summarize 
the equations to simplify their implementation. 

The volume element and the Laplacian are 


dV =a 3 (£ 2 — rj 2 )d^dr]d(j), 
1 


(Al) 


A = 


a 2 {e~V 2 ) 


V-i4 + £(i-V 
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e-v 2 
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The action of T = — A/2 on 'F m becomes 
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(A3) 
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(A7) 


The matrix elements of the dipole operator in ^-direction 
are 




ttfimm' $ij • 


(AS) 


The term stemming from the volume element, Eq. m, 
is canceled by the normalization factor of the basis func¬ 
tions, Eq. ra. 


2. Kinetic energy 


The factor (£ 2 — r) 2 ) 1 [e.g. Eq. (IA3I) ] is canceled by the 
volume element, which avoids numerical problems due to 


The matrix elements for the kinetic energy are evalu¬ 
ated using integration by parts and the DVR quadrature 
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3. 


Interaction energy 
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(A10) 


The Coulomb interaction of the electrons in Eq. © can 
be decomposed using the well-known Neumann series in 
prolate spheroidal coordinates [90, [91] : 


1 
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(A12) 


with £< = min(£i,£ 2 ) and £> = max((i,( 2 ). are 
the spherical harmonics and P z m (QJ 71 ) are the (irreg¬ 
ular) associated Legendre functions [78I . [92]. For their 
computation for £ > 1, the code from Ref. |93| was used. 

The 77 - and (^-dependent parts are evaluated straight¬ 
forwardly using the properties of the basis functions. Be¬ 


cause the Legendre functions exhibit a singularity for 
£ —> 1, the ^-dependent parts are not evaluated by DVR 
quadrature but by solving the corresponding differen¬ 
tial equation of the Green’s function Pj M \^<)Q\ M ^ (^>) 
[671 . l 68 l |. The final expression for the integrals is then 
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(A14) 


is the value of the last (excluded) FEDVR grid point, 
and Ui are the quadrature weights of the FEDVR grid 
points. In the implementation, the Legendre polynomials 
in 77 and the term in brackets are precomputed and stored 
in arrays so that the actual computation of the matrix 
elements is just summing up the product of three array 
values. The matrix elements are symmetric in £ and 77 , 
which can be exploited as well. 

Since the integration over 77 is still handled by usual 
quadrature, the maximum used value of l in the series 


expansion, Z max , should not be too large such that the 
integration kernel for the 77 -dependent part is not a poly¬ 
nomial of degree 2N V — 1 any more (N v is the number 
of functions in 77 ). Numerically, however, the results are 
not very sensitive to the choice of Z max and setting it to 
the number of used basis functions in 77 leads to good 
results MM- 
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Appendix B: Integral transformation 


The usage of a DVR-based basis with its diagonal¬ 
ly of potential and, to some extent, interaction ma¬ 
trix elements and the utilization of a partially rotated 
basis allows for a massive reduction of the usual scal¬ 
ing relationships in quantum-chemical algorithms. In 
the following, the most crucial ones, namely the gen¬ 
eration of the Fock matrix in Hartree-Fock and the in¬ 
tegral transformation to the basis of molecular orbitals 
are shown. For convenience, we use the “chemist’s” nota¬ 
tion [52] of the electron-electron repulsion integrals; i. e. , 
{ik\r^\jl) = ( ij\kl ). 


1. Coulomb matrix in Hartree-Fock 


The bottleneck in usual Hartree-Fock calculations is 
the generation of the Coulomb and the exchange matrix 
as parts of the Fock matrix. The elements of the Coulomb 
matrix are generated by 

Jjj = Du{ij\kl). (Bl) 

kl 

D is the density matrix, Eq. GZD- The scaling with the 
number of basis functions N b is O(N^) if no additional 
techniques like density fitting are used. For the exchange 
matrix, the summation runs over the integral ( ik\lj) and 
the optimization procedure is similar. 

Resolving the indices using Eqs. (TCH) and (IA13I) . 
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= E E Dkl fodidj ^a^ai &iiij fiiki-i 


mkikak mindi 


, 0 | rm-nijl 

A u mi-m j ,mi-mk i ‘‘i i ik,d i ,dk'‘ 


— A A \ ' \ ' n q| rm-rrijl 

— UtHdjViiij kl Hik,a.i,d k ’ 

m k =m B k i k a k 

-m max - (rrii - rrij), rrii- rrij < 0 
m.rTiax, else, 


l max) 


^maxj 

Wmax - ( rrii - rrij), else. 


rrii — rrij < 0 


(B2) 


(B3) 

(B4) 


The last two 5-symbols can be resolved in the overall 
loop creating the Coulomb matrix. Therefore, an overall 
scaling of less than is achieved for the construction 

of J, where is the number of basis functions in </>. The 
diagonalization of the Fock matrix is then the bottleneck 
of the SCF procedure. 


2. Integral transformation to the partially rotated 

basis 


The structure of the coefficient matrix for a partially 
rotated basis is 


c = 




(B5) 


where C rot is a dense matrix of size N IOt x N rot for the 
rotated block of basis functions. 


a. Integrals between rotated orbitals 

In the following, orbitals indexed with p , q, r or s de¬ 
note rotated orbitals (from the inner spatial region) and 
those indexed with a, b, c or d denote nonrotated orbitals 
from the inner region that have to be rotated. The 
rotated orbitals are unitarily transformed by the real¬ 
valued coefficient matrix C, Eq. (IB5I) . The integrals in 
the rotated molecular orbital frame are then |[52| 

N Iot 

(p\h\q) = Ca P C bq (a\h\b ), (B6) 

ab 

N rot 

C pq\rs ) = C ap C bq C cr C ds (ab\cd). (B7) 

abed 

Because of the structure of the coefficient matrix, 
Eq. (IB5I) . the sum runs only over all rotated orbitals, 
if { p , q, r, s} are all themselves rotated orbitals. It is well 
known that the formally 0(7V r 8 ot )-scaling transformation 
can be massively reduced by employing partial transfor¬ 
mations ElMieSl, i. e., first transforming and storing 
the fourth orbital to (ab\cs), then, transforming the third 
orbital to (a6|rs) and so on. This scales as 0(N^ ot ) but 
twice the memory is needed for storing the intermediate 
transformed integrals. Doing only two two-index trans¬ 
formation steps saves a considerable amount of memory 
but the scaling becomes worse. However, this is some¬ 
times favorable A similar procedure can be applied 
for the one-electron integrals, Eq. (IB6I) . 

In our case however, the transformation can again be 
sped up massively, as exemplified by the calculation of 
the Coulomb matrix, see Eq. (IB2I) . Because of the diag¬ 
onally in £ and (f>, it is beneficial to first transform the 
last two indices like in Eq. (IB2I) . where D^i has to be 
replaced with Ck a Ci b . Note that, although the integrals 
are real-valued, the orbitals for TO max ^ 0 are not: 

I rm,i 1-rrij , rm k « \ 

\J iidi J ijdj \ J ika^J iiai ) 

(ij\lk) 

/ rrrii r-rrij i rm k mi \ 

\J iidi J ijdj \J ikCLk ha-i ' 

Hence, only the following symmetries hold: 

(ij\kl) = (kl\ij) = (ji\lk) = ( lk\ji ) ± (lk\ij) ^ ... 

(BH) 


= (ij\kl) ^ (ij\lk) (B8) 

= ( I f™k \ /OQ\ 

iidiJ ijdjt’f iiaiJ \ u ) 

= ( f I fr m * fr m M (BIO) 

W lidi J ijdj \J iiai J ik a k ' ' / 
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Because the two-index transformed integrals have even 
less symmetry, _/V r 3 ot x (2m max + 1) elements need to be 
stored for them. The transformation of the remaining 
two indices are done using partial transformations, as 
usual in quantum chemistry (see above). For large non- 
rotated bases, the fully transformed integral tensor re¬ 
quires too much memory and the latter transformation 
is done on the fly: 

-Wrot ^max 

(j>q\rs) = T, Y, C ap C bq (ab\rs). (B12) 

d—\ 7Tlb = Hl m a X 

The index b on the right-hand-side is constructed by m b 
and the indices for the functions in f and 77 of the index 

a. The product C ap C bq can be precalculated which sim¬ 
plifies and accelerates the summation. Transforming the 
last two indices on the fly requires N rot x (2m max + 1) 
operations, which is still better than the usual require¬ 
ment of iV r 1 2 ot operations for a non-DVR basis such that 
the overall scaling for the complete two index transfor¬ 
mation is then lV r 5 ot x (2m max + 1) instead of the best 
achievable iV r 4 ot scaling. 

b. Integrals between nonrotated orbitals and rotated orbitals 

If all basis functions are nonrotated DVR functions, 
the formulas for the one- and two-electron integrals, 
Eq. (IA13I) . can be applied directly. The crucial point 
is the efficient implementation of the integrals between 
both nonrotated and rotated basis functions. A lot of 
simplifications come from the structure of the coefficient 
matrix, Eq. (IB5I) . which is diagonal if both basis func¬ 
tions are nonrotated functions and zero if one function 
is a rotated and another a nonrotated function. In the 
following, nonrotated basis functions are underlined. 

The mixed one-electron integrals are then: 

N b 

(p\h\q) = (q\h\p) = y C ap C bq (a\h\b), (B13) 

ab 

C bq = S b q, (B14) 

Mot 

=> (p|%> = y C ap {a\h\q). (B15) 

a 

Since the potential is diagonal, (p|V|< 7 ) = 0. The struc¬ 
ture of the kinetic energy matrix [diagonality for <f> func¬ 
tions, see Eq. (IA9I) . and a banded sparsity pattern] can 
be exploited as well. 
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